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GEOMETRICALLY NONLINEAR ANALYSIS OF 
LAMINATED ELASTIC STRUCTURES 

J. N. Reddy 

Virginia Polytechnic Institute and Stae University 
Blacksburg, VA 24061 

(Final Technical Report on the NASA Grant NAG-3-208) 

SUMMARY 

The research performed under the present grant deals with the 
analysis of laminated composite plates and shells that can beused to 
model automobile bodies, aircraft wings and fuselages, and pressure 
vessels among many others. The finite element method, a numerical 
technique for engineering analysis of structures, Is used to model the 
gometry and approximate the solution (l.e., displacements, stres SSS uud 
natural frequencies). Various alternative formulations for analyzing 
laminated plates and shells are developed and their finite element 
models are :ested for accuracy and economy in computation. These 
Include, the shear deformation laminate theory and degenerated 3-D 

elasticity theory for laminates. The present results obtained for stale 

\ 

transient and natural vibration olf laminated plates and shells are very 
accurate when compared to existing theories (e.g., classical lamination 
theory). Many of the results obtained during this Investigation should 
serve as references for future investigations by designers and 
experimental ists. 

Preliminary investigations were also initiated during this research 
program In two other related areas: the development of a refined shear 

deformation theory, and nonlinear constitutive models for composite 
plates. These Investigations are currently progressing under other 
sponsorships. 
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INTRODUCTION 

This final technical report contains three parts: Part 1 deals 

with the 2-D shell theory and Its finite element formulation and 
applications. Part 2 deals with the 3-D degenerated element. These two 
parts constitute the two major tasks that were completed under the 
grant. An other related topic that was initiated during the present 
Investigation s the development of a nonlinear material model. This 
topic Is briefly discussed in Part 3. To make each part self-contained, 
conclusions and references are Included in each part. In the interest 
of brevity, the discussions presented here are relatively brief. For 
details and additional topics see the journal articles 1 and 10 for Part 
1, articles 2 and 6 for Part 2 and article 11 for Part 3. 
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P ART 1 

GEOMETRICALLY NONLINEAR ANALYSIS OF LAMINATED SHELLS 
INCLUDING TRANSVERSE SHEAR STRAINS 


J. N. Reddy and K. Chandrashekhara 

(A cond&Mzd v&xA.ion ofa tka> papeA u to appear in AJ AA JouAnal . 19841 


SUMMARY 

The paper contains a description of a doubly curved shell finite element 
for geometri cal ly nonlinear (in the von Karman sense) analysis of laminated 
(doubly-curved) composite shells. The element is based on an extension of the 
Sanders shell theory and accounts for the von Karman strains and transverse 
shear strains. The numerical accuracy and convergence characteristics of the 
element are further evaluated by comparing the present results for the bending 
of isot opic and orthotropic plates and shells with those available in the 
literature. The many numerical results presented here for tne geomerti ca 1 ly 
nonlinear analysis of laminated composite shells should serve as reference for 
future investigations. 

INTRODUCTION 

Laminated shells are finding increased application in aerospace, automo- 
bile and petrochemical industries. This is primarily due to the high stiff- 
ness to weight ratio, high strength to weight ratio, and less machining and 
maintenance costs associated with composite structures. However, the analysis 
of composite structures is more complicated when compared to metallic struc- 
tures, because laminated composite structures are anisotropic and character- 
ized by bending-stretching coupling. Further, the classical shell theories, 
which are based on the Kirchhoff-Love kinematic hypothesis (see Naghdi [1] and 
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Bert [2]) are known to yield deflections and stresses in laminated shells that 
are as much as 30% in error. This error is due to the neglect of transverse 
shear strains in the classical shell tneories. 

Refinements of the classical shell theories (e.g., Love's first approxi- 
mation theory [?]) for shells to include transverse shear deformation have 
been presented by Reissner Sanders [7] presented modified first- and 

second-approximation theories that removed an inconsistency ( non\ani shi ng of a 
small riuid-body rotations of the shell) existed in Love's first-approximation 
theory . 

The first thin shell theory of laminated orthotropic composite shells is 
due to Ambartsumyan [3,9]. In L.iese works Ambartsumyan assumed that the indi- 
vidual orthotropic layers were oriented such chat the principal axes of mate- 
rial symmetry coincided with the principal coordinates of the shell reference 
surface. Dony, Pister, and Taylor [10] presented an extension of Donnell's 
shallow shell theory [11] to thin laminated shells. Using the asymptotic in- 
tegration of the elasticity equations, Widera and Chung [12] derived a first- 
approximation theory for the unsymmetric deformation of nonhomoyeneous , aniso- 
tropic, cylindrical shells. This theory, when specialized to isotropic mate- 
rials, reduces to Donnell's shell theory. 

Tne effects of transverse shear deformation and thermal expansion through 
the snell thickness were considered by Zukas and Vinson [13]. Dong and Tso 
[14] constructed a laminated orthotropic shell theory that includes transverse 
shear deformati on . This theory can be regarded as an extension of Love's 
first-approximation theory [3] for homogeneous isotropic shells. Other re- 
fined theories, specialized to anisotropic cylindrical shells, were presented 
by Whitney and Sun [15]. and Widera and Logan [16,17]. 
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The finite-element analysis of layered anisotropic shells, all of which 
are concerned with bending, stability, or vibration of shells, can be found in 
the work? of Schmit and Monforton [18], Panda and Natarajan [19], Shivakjmar 
and Krishna Murty [20], Rao [21], Siede and Chang [22], Hsu, Reddy, and Bert 
[23], Reddy [24], and Venkatesh and Rao [25,26], Recently, Reddy [27] extend- 
ed the Sanders theory to account for the transverse shear strains, and pre- 
sented exact solutions for simply supported cross-ply laminated shells. All 
of these studies are limited small displacement theories and static analyses. 

In the present paper, an extension of the Sanders shell theory that ac- 
counts for the shear deformation and the von Karman strains in laminated an- 
isotropic shells is used to develop a displacement finite element model for 
the bending analysis of laminated composite shells. The accuracy of the ele- 
ment is evaluated by comparing the results obtained in the present study for 
isotropic and orthotroDic plate and shell problems with those available in the 
literature. Numerical results for bending analysis of cylindrical and doubly- 
curved shells are presented, showing the effect of radius-to-thickness ratio, 
loading, and boundary conditions on the deflections and stresses. 

A REVIEW OF THE GOVERNING EQUATIONS 

Consider a laminated shell constructed of a finite number of uniform- 
thickness orthotropic layers, oriented arbitrarily with respect to the shell 
coordinates U lt £otC)* ^he ortho 9 ona l curvilinear coordinate system 
Ui*S2>C) chosen SuCf1 the Z>i~ anc * curves are lines of curvature on 

the midsurface £=Q, and ^-curves are straight lines perpendicular to the sur- 
face c~ 0 (see Fig. 1). A line element of the shell is given by (see Reddy 
[27] 


(dS) 2 = [[ i + C/R^a^J 2 + [(1 + C/R 2 )a 2 d^2] 2 + (dc) 2 


( 1 ) 
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where a i and Rj (1 » 1,2) are the surface metrics and radii of curvature, 
respectively. In general, and R n - are functions of ^ . only. For the doubly 
curved shells considered in the present theory, a 1 and R^ are constant. 

The strain-displ acement equations of the shear deformable theory of 
doubly-curved shells are given by 


= e 1 + 

- 0 ± r 

E £ “ E^ £ X £ 


where 


-0 s 


e 5 s e 5 
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Here u^ denotes the displacements of the reference surface along = z) 

axes, and are the rotations of the transverse normals to the reference 
surface. In Love's first-approximation theories the parameter c Q is taken to 
be zero, and it is introduced only in the Sanders theory. 
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The stress-strain relations, transformed to the shell coordinates, are of 
the form 

(o} = [Q] {e} (4) 


(k) 

where Qi . ' are the material properties of k-th layer. 

* J 

The principle of virtual work for the present problem is given by 


n - V r i 00* 4. ( k ) . . (k). . <k). . (k). 

0 ~ L J / + a 2 6e 2 °6 6e 6 + a 4 6e 4 % 6e £ 


k=l q 


k-1 


- qfiu^jciia^di^d^dC 


5) 


= / +• + ^g e 5 + + ‘^2 I ^ I? 2 + ^6^ K 5 + ^l^ e 5 

Q 

+ Q 2 6e3 - qfii^Ja^d^d!^ (6) 


where q is the distributed transverse load, and are the stress and 
moment resultants, and is the shear force resultant: 

L c k 

(N. ,M ) = l f a.{l,C)dC , i = 1,2,6 
k=1 <k-l 
L C k 

Qi - l K? / a,dc , i = 4,5, (7) 

k=1 «k-l 

where (i = 1,2) are the shear correction factors (taken to be = 

5/6), and (c k ^C^) are the ^-coordi nates of the k-th layer, and l is the 
total number of layers in the laminated shell. 

It is informative to note that the equations of equilibrium can be 
derived from Eq. (6) by integrating the displacement gradients in e° by parts 
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and setting the coefficients of 6u^ (i = 1,2,3) and 6^ ( i =1 ,2 ) to zero 
separately. We obtain [with c Q = j (|~ - ^-) and dx^ = a^d^] 


a Qi 

5x^ + axj ' N 5 ■ c oV + “ 0 


8N, Q 2 

^ < N 6 + C °V * 8^ + Rj * 0 


SQl 5Q 2 N N 2 


ax^ ax 2 


- (<r + T- 0 + W K) = 0 


■ R i r 2 


dM aM 

4 + 4 ' Ql = ° 


5M fi &M- 

8** + axf ' Q 2 " 0 


( 8 ) 


where 


au_ au 

’af^l Fxf* N 6 57- 


*(“3) ■ 57T (" 


> N 
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2 ax. 


(9) 


The resultants ( N ^ , Q i ) are related to ( e °, ) (i,j = 1,2,6) by 


N i ■ VI + B ij-j 

M i ■ V° * D ij-j 

^2 = A 44 £ 4 + A 45 e 5 

Q 1 = A 45 e 4 + A 55 g 5 


( 10 ) 

(ID 


Here Ay , By ard Dy ( i , j = 1,2,6) denote the extensional, flexural' 
extensional coupling, and flexural stiffnesses of the laminate: 


L 


(A^.B.-.D.J = l J Q^ } (l,C,C 2 )dC (U = 1,2,6) 


1J* 1J 1J 


k=l C 


k-l 


ij 


( 12 ) 
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( A 44 ,A 45’ A 55^ = ,K l K 2 Q 45^ K 2 q 55^ dC 

The boundary conditions, derived from the virtual work statement, involve 
specifying either the essential boundary conditions ( EBC) or the natural 
boundary conditions (NBC): 


EBC NBC 
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(13) 


where (n^,n 2 ) denote the direction cosines of the unit normal on the boundary 
of the midsurface of the shell. 

The exact form of the spatial variation of the solution of Eqs, (8)-(13), 
for the smal 1 -displacement theory, can be obtained under the following condi- 
tions (see Reddy [27]): 

(i) Symmetric or antisymmetric cross-ply laminates: i.e., laminates 

with 

A 16 = A 26 = B 16 = B 26 = °16 = °26 “ A 45 = °* 

(ii) Freely supported boundary conditions: 

N 1 (0,x 2 ) = Nj(a,x 2 ) = M 1 (0 ,x 2 ) = M 1 (a,x 2 ) = 0 
u 3 (0,x 2 ) = u 3 (a,x 2 ) = u 2^ 0,x 2^ 3 u 2^ a,X 2^ " 0 
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N 2 (x 1 .0) = N 2 (x r b) = M 2 (x r O) = M 2 (x r b) = 0 

u 3 (x 1 ,0) - u 3 (x 1? b) = u^Xj.O) = UjUj.b) = 0 

<J) 2 ( 0,x 2 ) = t)> 2 ( a ,x 2 ) => ^(x^O) = ^(x r b) = 0 (15) 

where a and b are the dimensions of the shell middle surface along 

the Xj and x 2 axes, respectively. The time variation of the load 
does not influence the spatial form of the solution. 

Note that the exact solution can be obtained only for cross-ply laminated 
shells with simply supported boundary conditions. For general lamination 
schemes, exact solutions are not available to date. 


F INITE-ELEMENT MODEL 

A typical finite element is a doubly-curved shell element in the x^x 2 - 

( e ) 

surface. Over the typical shell element g v , the displacements 
(u^ »u 2 > u 3 ,<t>^ , 4 > 2 ) are interpolated by expressions of the form. 


u 1 = .1 ufyjUj ,x 2 ) , i = 1 ,2,3 


j"l 


« i 

$,= I »x„) , i = 1 .2 

1 3=1 2 


(16) 


where 4 . . are the interpolation functions, and ulj and are the nodal values 
of u j and $ respectively. For a linear isoparametric element (N = 4) this 
interpolation results in a stiffness matrix of order 20 by 20. For a nine- 
node quadratic element the element stiffness matrix is of order 45 by 45. 

Substitution of Eq. (21) into the virtual work principle, Eq. (9) yields 
an element equation of the form 
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[K(a)] (A) = {F} (17) 

where {a} = {{uj}, {u 2 } , {u 3 } , » {4>^ } } T » [K] the elemert stiffness matrix, 

and { F} is the force vector. In the interest of brevity, the coefficients of 
the stiffness matrices are included in Appendix I, 

The element equations (17) can be assembled, boundary conditions can be 
imposed, and the resulting equations can be solved at each load step. Note 
that the stiffness matrix [K] is a function of the unknown solution vector 
(a); therefore, an iterative solution procedure is required for each load 
step. In the present study, we used the direct iteration technique, which can 
be expressed as 

[K(w r )iw r+l - in us) 

p 

where {a} denotes the solution vector obtained in the r-th iteration (at any 
given load step). At the beginning of the first load step, we assume that 
(A}° = {0} ano obtain the linear solution at the end of the first iteration. 
The solution ottained at the end of the r-th iteration is used to compute the 
stiffness matrix for the (r+l)-th iteration. At the end of each iteration 
(for any load step), the solutions obtained in two consecutive iterations are 
compared to see if they are close enough to terminate the iteration and to 
move on to the next load step. The following convergence criterion is used in 
the present study: 

[ J Uj-4 +1 | 2 / I |a?I 2 ] 1/z < 0.01 (19) 

i=l 1 i=l 1 

where N is the total number of unknown generalized displacements in the finite 


element mesh. 
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To accelerate the convergence, a weighted average of the solution from 
last two iterations are used to compute the stiffness matrix: 

[K( Y {A} r-1 + (1 - Y ){A} r )]{A} r+1 = {F} (20) 

where y is the acceleration parameter, 0 < y < 1. In the present study a val 
ue of 0.25 - 0.35 was used. 


NUMERICAL RESULTS 

Here we present numerical results for Lome sample problems. To illus- 
trate the accuracy of the present element, first few examples are taken from 
the literature on isotropic and orthotropic shells. Then results (i.., de- 
flections and strsses) for several laminated shell problems are presented. 

The results for laminated shells should serve as references for future inves- 
tigations . 

All of the results reported here were obtained using the double-precision 
arithmatic on an IBM 3081 processor. Most of the sample problems were an- 
alyzed using a 2 x 2 uniform mesh of the nine-node (quadratic) i soparametric 
rectangular element. 

1 . 8ending of a simply supported plate strip (or, equivalently, a beam) under 
uniformly distributed load. 

The problem is mathematically one-dimensional and a,i analytical solution 
of the problem, based on the classical theory, can be found in Timoshenko and 
Womowsky-Krieger [28]. The plate length along the y-coordinate is assumed to 
be large compared to the width, and it is simply supported on edges parallel 


to the y-axis. The following simply supported boundary condition. are 
assumed : 

3 $2 = 0 along edges x = - 127mm (21) 

All inplane displacement degrees of freedom are restrained. A 5 x 1 mesh of 
four-node rectangular elements in the half plate is used to analyze the prob- 
lem. The data and results are presented in Fig. 2. The present result is in 
good agreement with the analytical solution. 

2 . Clamped square plate under uniform load. 

Due to the biaxial symmetry, only one quadrant of the plate is modelled 
with the 2x2 mesh of nine-node elements (4x4 mesh of linear elements give 
almost the same result). Pertinent data and results are presented in Fig. 3 
for side to thickness ratios a/h = 10 and 500. The result for a/h = 500 is in 
agreement with the results of Way [29]. The difference is attributed to the 
fact that the present model includes the inplane disp'acement degrees of free- 
dom and transverse shear deformation. 

Figure 4 contains transverse deflection versus load for clamped ortho- 
tropic, cross-ply, and angle-ply plates. The lamina properties are 

E 1 = 25 x 10 4 N/mm 2 , E 2 » 2 x 10 4 N/mm 2 , G 12 = G 13 = 10 4 N/mm 2 

G 23 3 0.4 x 1G 4 N/mm 2 , Vj 2 = 0.25. 

For the same total thickness the clamped orthotropic square plate is stiffer 
than both two-layer angle-ply and cross-ply plates. 









Load 

% 

t lO'^N/mm 2 ) 



Deflection, ~w (in mm) 


Figure 4. Bending of clamped orthotropic and laminated 
square plates under uniform load. 
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3 . Simply supported, isotropic spherical shell under point load. 

The pertinent data of the 3 nel 1 is shown in Fig. 5. A uniform mesh of 2 

by 2 quadratic elements is used in a quadrant. The effect of three types of 

simply supported conditions on the center deflection and center normal stress 
is investigated: 

SS-1: u = w = 4^=0 at y = b; v = w = i> 2 =0 at x = a 

SS-2 : u = v = w = 4> t =0 at y = b; u = v = w = <p 0 = 0 at x = a (22) 

SS-3: v = w = $ =* 0 at y = b; u = w = = 0 at x = a 

Table 1 contains the results for the three boundary conditions. It is clear 
from the results that all three boundary conditions give virtually the same 
results for a/h - 150, and d fer significantly (especially SS-1 differs from 
both SS-2 and SS-3) for a/h = 16. Thus, the effect is more in thick shells 
than in thin shells. The stress o x shown in Fig. 5 is evaluated a* point x = 
y = 1.691" in the top layer 

4. Simply support e d Isotropic cylindrical shell under point load. 

The geometry and finite-element mesh of the shell are shown In Fig. 6. 

Once again, the effect of various simply supported boundary conditions (22) on 
the deflections and stresses for the problem is investigated using a uniform 
mesh of 2 x 2 quadratic elements. The results are presented in Table 2. For 
the geometry and loading used here (R = 2540, a = 254, h = 12.7), 'he boundary 
conditions have very significant effect on the solution. Boundary conditions 
SS-2 and SS-3 give almost the same results whereas SS-1 gives about 2-1/2 
times the deflection given by SS-2 or SS-3 boundary conditions. 
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Table 1. Effect of various simply supported boundary conditions on the center 
deflections and normal stress in spherical shells und^r point load 
( E = 10 7 psi , v = 0.3) . 


Load 

P/h 2 

Solution 

a/h=16Q SS 

-1 

a/h=16 

SS- 

a/h=160 

■2 

a/h=16 

SS- 

a/h=l60 

■3 

a/h=16 

4,000 

-w* 

0.0155 

_ 

0.0152 


0.0152 



-a* 

893 

- 

984 

- 

894 

- 

8,000 

-w x 

0.0329 

0.0349 

0.0324 

0.0255 

0.0324 

0.0258 


_CJ y 

1,880 

6,535 

1 ,882 

6,015 

1,882 

6,031 

12,000 

X 

-w 

0.0529 

- 

0.0522 

- 

0.0521 

- 



2,980 

- 

2,985 

- 

2,986 

- 

16,000 

-w 

0.0760 

0.0793 

0.0752 

0.0520 

0.0751 

0.0525 


" a y 

4,220 

13,230 

4.228 

12,200 

4,229 

12,240 

20,000 

-w x 

0.1038 


0.1028 

- 

0.1027 

- 



5,657 

* 

5,671 

- 

5,672 

- 

24 ,000 

-w 

0.1364 

0.1083 

0.1354 

0,0792 

0.1353 

0.0800 



7,268 

20,110 

7,289 

18,500 

7,291 

18,550 

28,000 

X 

-w 

0.1761 

~ 

0.1752 

- 

0.1751 

- 


”0 

9,128 

- 

9,160 

- 

9,162 

- 

32,000 

-w x 

r\ hdh a 
U 

0.1472 

U o £ CCf 

r\ i ,■> 1 n 

U *1U/ c 

0.2227 

0.1083 


•°x 

11,180 

27,170 

11,220 

24,930 

11,230 

25,000 

* w( 0,0) 

» <? X (A,A); 

A = 1.691 







Table 2. 

Effect of various types of simply supported boundary conditions on 
the deflections and stresses of anisotropic cylindrical shell jnder 
point load. 

Load ,P 

SS- 

1 

SS-2 


SS- 

3 

(N) 

-w ( mm ) 

-Oy ( N/mm 2 ) 

-w 

" a y 

-w 

~ a y 

250 

2.5804(2) 

2.868 

0.6544(4) 

1.706 

0.6698(4) 

1.706 

500 

5.1626(2) 

5.713 

1.3533(4) 

3.478 

1.3843(4) 

3.477 

750 

7.7343(2) 

8.506 

2.1057(4) 

5.327 

2.1522(4) 

5.321 

1,000 

10.278(2) 

11.210 

2,9234(4) 

7.265 

2.9855(4) 

7.242 

1,250 

12.733(2) 

13.80 

3.8241(4) 

9.312 

3.9017(4) 

9.288 

1.500 

15.204(2) 

16.25 

4.8349(4) 

11.50 

4.9279(4) 

11 .46 

1,750 

17.560(2) 

18.560 

6,0331(5) 

13.91 

6.3 423(5) 

13.85 

2,000 

19.843(2) 

20.730 

7.5316(6) 

16.66 

7.5610(6) 

16.57 



Deflection 
w( in. ) 



0 80 160 240 320 

load, P (lbs) 


Figure 5. Bending of a simply supported (SS-3), 
isotropic, spherical shell under point 
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5 . Clamped Isotropic cylindrica shell under uniform loading. 

Figure 7 contains the pertinent data and results for a clamped cylindri- 
cal shell (isotropic) subjected to uniform load. The results are compared 
with those obtained by Dhatt [30], 7 he agreement is very good. 

6. Clamped orthotropic cylindrical shell subjected to internal pressure. 

Figure 8 contains the geometry and plots of center deflection and center 
stress versus the internal pressure for the problem. The orthotropic material 
properties used in the present study are: 

3 7,5 * 10 6 psi , E ? = 2 x 10 6 psi , G^ = = G 0 o = 1 .25 x 10 6 p si 

\>12 = 0,25 (23) 

The ^resent result, obtained using the 2x2 mesh of quadratic elements, is in 
excellent agreement with that obtained by Chang and Sawamiphakdi [31]. 

7 . Nine-layer [0 o /90 o /0° . . ,/Q°] cross-ply spherical shell subjected to 
uniformly distributed load, 

Th<- following geometrical data is used in the analysis (with SS-3 boundary 
conditii ' ■ ) : 

= R 2 = R = 1,000 in., a = b = 100 in., h = 1 in. (24) 

Individual layers are assumed to be of equal thickness (h^ = h/9), with the 
zero-degree layers being the inner ana outer layers. The following two sets 
of orthotropic-material constants, typical high modulus graphite epoxy materi- 
al (the ratios are more pertinent here), for individual layers are used: 
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Mat o-l : E^ = 25 x 10 6 psi, 3 10 6 psi , G 12 = s 0*5 x 10 6 psi 

G 23 = 0.2 x 10 6 psi, v 12 = 0 .25 (25) 

Mat, -2: Ej = 40 x 10 6 psi, E 2 = 10 6 psi, G 12 3 G 13 = 0.6 x 10 6 psi 

G 23 = 0.5 x 10 6 psi, v 12 = 0.25 (26) 

Figure 9 contains plots of center deflector (w/h) versus the load parameter 

(P = q Q R 2 /E 2 h 2 ) for the two materials. Shell constructed of Material 1 
deflects more, for a given load, than the shell laminated of Material 2 
(because Material 2 is stiffer), and consequently experiences greater degree 
of nonlinearity. Note that the difference between the nonlinear deflections 
of the two shells increase nonlinearly, indicating that the shell made of 
Material 2 can take much more (ultimate) toad than apparent from the ratio of 
moduli of the two materials, E^/E^. 

8 . Effect of various simply-supported boundary conditions on the deflections 
of two-layer cross-ply spherical shells under uniform load. 

As pointed out in Problems 3 and 4, the transverse deflection is sensi- 
tive to the boundary conditions on the inplane displacements of simply sup- 
ported shells. To further illustrate this effect for laminated shells, a set 
of four types of boundary conditions are used, and the results are presented 
in Table 3. Here SS-4 has the following meaning: 

w = ^ = 0 on x = a 
w=<f> 2 =0ony = b 


SS-4 


(27) 



Deflection, ( -w/h ) 


3. Bending of nine-layer cross-ply 
[0°/90°/0°/. . . ] spherical shell 
subjected to uniformly distributed 
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Table 3. Effect of various simply-supported boundary conditions on the trans- 
verse deflections of cross-ply [0°/90°] spherical shells under 
uniform load (Material 1; shell dimensions are the same as those in Fig. 10) 


% 

(psi) 


-w 

(In.) 


SS-1 

SS-2 

SS-3 

SS-4 

0.50 

0.3344 

0.04246 

0.04257 

0.4592 

0.75 

0.5757 

0.06599 

0.06617 

0.8255 

1.00 

0 .9485 

0.09144 

0.09171 

1.3845 

1.25 

1.6529 

0.11926 

0.11966 

1.9589 

1.50 

2.2826 

0.15008 

0.15063 

2.3597 

1.75 

2.6421 

0.18478 

0.18556 

2.5951 

2.00 

2.8499 

0.22473 

0.22584 

2.8074 

2.25 

3.0764 

0.27425 

0.27593 

3.0284 

2.50 

3.2432 

0.33534 

0.33795 

3.1948 

2.75 

3.4214 

0.42970 

0.43487 

3.3719 
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Once again we note that SS-2 and SS-3 give almost the same deflections. 
Boundary conditions SS-1 and SS-4 give deflections an order of magnitude high- 
er than those given by SS-2 and SS-3. Thus, boundary conditions SS-2 and SS-3 
make the shell quite stiffen. 

9. Two-layer cross-ply [Q°/90°] and angle-ply [-45°/45°], simply-supported 
(SS-3) spherical snells. 

Figure 10 contains the pertinent data and results (with different scales) 
for the cross-ply and angle-ply shells (of Material 2). It is interesting to 
note that the type of nonlinearity exhibited by the two shells is quite dif- 
ferent; the cross-ply shell gets softer whereas the angle-ply shell gets 
stiffer with an increase in the applied load. While both shells have bending- 
stretching coupling due to the lamination scheme = " 8 n nonzero for the 

cross-ply shell and B-^ and E i ^ are nonzero for the angle-ply shell), the 
angie-ply experiences shear coupling that stiffens the spherical shell rela- 
tively more than the normal coupling (note that, in general, shells get softer 
under externally applied inward loau). 

Figure 11 contains plots of center deflection, normal stress (-a ) and 
shear stress (cr ) at x = y = 5,233" versus load for two-layer cross-ply 
(0°/90°) spherical shell (Material 1) under point load at the center of the 
shell. The nonlinearity exhibited by the stresses (especially a ) is less 

yz 

compared to that exhibited by the transverse deflection. 

10. Two-layer clamped cylindrical shells under uniform loads. 

Figures 12 and 13 contain results (i.e., w, a , a versus load) for 

y xz 

cross-ply [0°/90°] and angle-ply [-45°/45°j clamped cylindrical shells under 
uniform load. The load-deflection curve for the cros-ply shell resembles that 


Deflection, -w (10“ in) -* — angle-ply 



Deflection, w (in in. ) 


Figure 10. Sending of two-layer cross-ply and angle-ply, 
simply supported (SS-3) spherical shells under 
uniform load. 



P( 10 3 1 bs ) 



Deflection and Stresses, -w,-a ,-o 

y yz 


Figure 11. Bending of a cross-ply [0 7 90°] spherical shell 

(SS-3, Material 1), under point load. (see Fig. 10 
for the shell dimensions) 
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20 

40 
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80 
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fiyure 13. Bending of a clamped cross-ply [0V90°] cylindrical 
shell under uniform load (Material 1) 
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of the isotropic shell in Fig. 7, but exhibits greater degree of nonlinearity 
(being stiffer). The angle-ply shell exhibits different type of nonlinearity 
(softening type) for all loads. 


11 . Quasi-isotropic, clamped, cylindrical shell under uniform load. 

Two types of quasi-isotropic clamped cylindrical shells are analyzed: 


Type 1: [0 J /45°/90°/-45°] sym> 

Type 2: C0°/+45 °/90 J $ym ^ 


(28) 


Material 1 properties are assumed for each lamina (8 layers). The geometric 
data and results are presented in Fig. 14. Compared to the results presented 
in Figs. 12 and 13, ten quasi-i sotropic shells have the 'near-inflection' 
point at higher loads; the load-deflection curve has essentially the same form 
as that of the cross-ply shell (see Fig. 12). 


CONCLUSIONS 

A shear-flexible finite element based on the shear deformation version of 
the Sanders' theory and the von Karman strains is developed, and its applica- 
tion to isotropic, orthotropic, and laminated (cross-ply and angle-ply) shells 
is illustrated via numerous sample problems. Many of the results, especially 
those of laminated shells, are not available in the literature and therefore 
should serve as references for future investigations. From the numerical com- 
putations it is observed that boundary conditions on the inplane displacements 
have significant effect on the shell deflections and stresses. Also, it is 
noted that the form of nonlinearities exhibited by different lamination schemes. 
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APPENDIX I 

Stiffness Coefficients : 

Let f -1^1 f = I ^3 
L6t T 1 " 2 2 ^ ' f 2 2 dx 2 

[Kll] = A n LS u ] + A 16 ([S12] + [S21]) + A 66 [S22] 

- C o( B I 6 (CS12] + CS21] ) + ZB 66 Csa23 ' C o° 66 [S22] ) + ^SOO] 

R 1 

[« 12 ] = A 12 [S12] + A 16 [S^] + A 26 CS22] + A 66 [S21] 

- S(V- s22] - c oB 66 [S»]J ♦ ^CSOO] 

[K“] = fjfAjjCS'l] + A 16 CS12] + [S21]) + A 6 «[S22]) 

+ f 2 (A 12 [S!2] + A 16 [Sl>] + A 26 [S22] t A 6g (S2‘]) 

* ^ * fl isCs 2 °3) * (A 12 [S‘»] * A 26 [S2«]) 

- c o(l?f Cs20] * CS20:I ) - 17 ('' 45 [S ° 2 J + 

- C o[ f l( B 16[ S21 ] * S 66 CS22]J ‘ f 2 (B 26 [S22] * 8 66 [S2>3)1 
[K‘A] . B u [S‘l] +B 16 ([S‘2] + Cs 21 j) + BggES 22 ] 
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[K 35 ] = 8 12 [S‘2] * B, 6 [Sl‘j *■ B^CS 22 ] ' B 66 [S 2 ‘] 

- W S22] + °66 [S21] ) ' ^ V S °° ] 

[K 23 ] = [K 32 ] T 

[K 22 ] = A 22 [S 22 ] + A 26 ([S‘ 2 : * [S 21 ] 3 * A m [S>1] + 2c 0 B 66 [S"] 

+ c o (8 26 ([S‘ 2 ] + [S 2 1]) + c o 0 66 [SU])-^-CS«»] 

R 2 

[K 23 ] = f^A^CS 23 ] + A, 6 [S 22 ] + A 16 [S«‘] *■ A 66 [S 32 ]) 

<• f 2 (A 22 [S 22 ] + A 26 (CS 21 ] ♦ lS 32 ] ) + A 66 [3Hj) 

+ ^ (A 12 i:s 2 "] * A 16 csi»]) + ^ (A 22 [S 2 «] * a 26 [s‘«]) 

+ c 0 * ^ csl0] - CA 44 CSO^ * A 45 [S» 3 ]) 

♦ St f i( B i6 [sll] + he^'l * V B 26 csl2] + 

[K 2 -] = b 12 [s 23 ] ♦ 3 26 [S 22 ] + B 16 [Sl>] + 3 SS [S12] 

+ c 0 (D 16 [SM] * °66 tSl2] ) - Rj A 45 CS °°] 

[K 25 ] . B 22 [S 22 ] + B 26 ([S 2 ‘] + [S 32 ]) ♦ B 66 [S 3 ‘] 

* c „( D 26 CsU] + D 66 [S ‘ 1] ) -r z V S °° ] 
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[K“] nl 
[K»] nl 
[k 33 i • 


[K 3 *] = 


[2K 13 ]jJ|l , NL = Nonlinear portion of the matrix 


= [2X23] 


NL 


a 45 Cs 12 ! + A„[S^] + A 14 [S-22] + A 45 [S21] 


55 


44 L 


+ 2ES“](A u fj + 2A 16 f,f 2 + A 66 f2 2 ) 

+ 2([S12] + [ S 21])[f2A 16 + (A 12 + A 66 )f 1 f 2 + f|A 26 ] 


+ 2[S 22 ](Aggf| + 2A 2g f 1 f 2 + A 22 f 2 ) 

. rcoo-iri- (hi + hi) + ±- (hl + hi)] 
+ [S°°][ R (-R- + R J + R l D + rJj 


'1 1 


'2 2 1 
A,, A, 




A - * A^~ A * 0 

([s o 2 ] ♦ 2[S 20 ])[fi( -l2. ♦ -|) + f 2 (-^- + Tr)] 


a 55 [s‘»] - a 45 [s 3 «] 


+ 2f 1 (B u [S il ] + B lg ([S 12 ] + [S 21 ]) + 3gg[S 22 ] ) 


+ 2f 2 (S 12 [S 2 n + B fifi [Sl 2 ] + B 2 g[S 22 ] + B lg CS“]) 


66 


B, , B. 


B,c B, 


t ( ^i. + ^3 + (-11 . -|l )[s o 2] 



[K 35 ] - A 45 csi»] * A 44 [S 3 «] 

* Bf^CsU] * b 66 [s 3 1] <■ a 16 csii] + B 26 [S 33 ]) 

+ 2f 2 (8 22 CS 33 ] * B 26 [S^] + B 26 [S 3 1] * B 66 [S“]) 

* ^)ts° 2 ] 

K 1 K 2 K 1 2 

['<41] „ [ K 14]T j [K 42] = [ K 24 ] T , [K43] nl = [K 24 ]^ 

[K 44 ] = D U [S“] + D 16 ([Sl*] + [S21]) + D 56 [S 22 ] ^ A 55 [S°°] 
[k 4 5] = D 12 [si2] + D lg CSH] + 0 26 [S 22 ] + D 66 [S21] + a 45 [s°°] 
[<51J = [K^f, [ K 52] = [K 25 ] T , [K53 ] nl = | [K35]J l 
[K54] = [k 4 5] T 

[k 55 ] = o 22 [s 22 ] + d 26 (Csi 2 ] + [s 2 !]) + o 66 [s“] + a 44 [s°°] 

^Linear " ^Unear 


where 
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It should be noted that although and f 2 are shown factored outside the 
matrices. In the evaluation of the coefficients by the Gauss quadrature fj and 
f 2 are considered as parts of the integrals. For example f 1 A 11 [S 11 ] is 
evaluated by 


/ f 1 A u 4 ) .4jdx 1 dx 2 
Q e 


, N N au- 

I l l 

1-1 J=1 u 3 1 


i^j 


»x 


W.W.detJ 
2 aZ J 1 0 


o 


where N is the number of Gauss points, Wj and Wj are the Gauss weights, Zj and 
Zj are the Gauss points, and J Q is the Jacobian of the transformation . 
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PART 2 


ANALYSIS OF LAMINATED COMPOSITE SHELLS 
USING A DEGENERATED 3-D ELEMENT 

W. C. Chao* and J. N. Reddy 
Department of Engineering Science and Mechanics 

(Tfu, i papoJi -U to appzcui in Int. JouAnal o ^ Num<zAtc.cil Mzth.ocU> in Engng . ) 

SUMMARY 

A special three-dimensional element based on the total Lagrangian 
description of the motion of a layered anisotropic composite medium Is 
developed, validated, and employed to analyze laminated anisotropic 
composite shells. The element contains the following features: 
geometric nonlinear I ty, dynamic (transient) behavior, and arbitrary 
lamination scheme and lamina properties. Numerical results of nonlinear 
bending, natural vibration, and transient response are presented to 
Illustrate the capabilities of the element. 

INTRODUCTION 

Composite materials and reinforced plastics are Increasingly used 
In automobiles, space vehicles, and pressure vessels. With the Increased 
use of fiber-reinforced composites as structural elements, studies 
involving the thermomechanical behavior of shell components made of 
composites are receiving considerable attention. Functional 
requirements and economic considerations of design have forced designers 
to use accurate but economical methods of determining stresses, natural 
frequencies, buckling loads etc. 


Graduate Research Assistant; presently at the University of Dayton 
Research Institute 
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Majority of the research papers In the open literature on shells 
Is concerned with bending, vibration, and buckling of J tropic 
shells. As composites materials are making their way Into many 
engineering structures, analyses of shells made of such materials 
becomes Important. The application of advanced fiber composites In jet 
engine fan or compressor blades and high performance aircraft require 
studies Involving transient response of composite shell structures to 
assess the capability of these materials under dynamic loads. 

Finite-element analysis of shell structures In the past have used 
one of the three types of elements: 1. a 2-D element based on a two- 

dimensional shell theory; 2. a 3-D element based on three-dimensional 
elasticity theory of shells; and 3. a 3-0 degenerated element derived 
from the 3-0 elasticity theory of shells. The 2-D shell theory Is 
derived form the three dimensional continuum field equations via 
simplifying assumptions. The simplifications require the Introduction 
of the static and kinematic resultants, which are used to describe the 
equations of motion. The unavailability of a convenient general 
nonlinear 2-D shell theory makes the 2-D shell element restrictive In 
Its use. The degree of geometric nonlinearity Included in the 2-D shell 
element Is that of the von Karman plate theory. In contrast to the 2-D 
shell theory, no specific shell theory Is employed In the 3-D 
degenerated element; Instead, the geometry and the displacement fields 
are directly discretized and Interpolated as In the analysis of 
continuum problems. 

Finite-element analyses of the large-displacement theory of solids 
are based on the principle of virtual work or the associated principle 
of stationary potential energy. Horrigmoe and Bergan [lj presented 
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classical variational principles for nonlienar problems by considering 
Incremental deformations of a contlnum. A survey of various principles 
In Incremental form Is presented by Wunderlich [2J. Stricklin et al. 

[3J presented a survey of various formulations and solution procedures 
for nonlinear static and dynamic structural analysis. The formulations 
include the pseudo force method, the total tagranglan method, the 
pdated Lagranglan method, ana the convected coordinate method. 

The only large-deflection analyses of laminated composite shells 
that can be found in the literature are the static analysis of Noor and 
Hartley [4] and Chang and Sawamlphakdl (51. Noor and Hartley employed 
the shallow shell theory with transverse shear strains and geometric 
nonl Inearltles tc develop triangular and quadrilateral finite 
elements. Chang and Sawamlphakdl presented a formulation of the 3-0 
degenerated element ter geometrically nonlinear bending analysis of 
laminated composite shells. The formulation Is based on the updated 
Lagranglan description and It does not Include any numerical results for 
laminated shells. 

From the review of the literature It Is clear that first, there 
does not exist any finite-element analysis of geometrically nonlinear 
transient response of laminated anisotropic shells, and second, the 3-D 
degenerated element Is not exploited for geometrically nonlinear 
analysis of laminated anisotropic shells. In view of these 
observations, the present study was undertaken to develop a finite- 
element analysis capability for the static and dynamic analysis of 
geometrically nonlinear theory of laminated anisotropic shells. A 3-0 
degenerated element with total Lagranglan description Is developed and 
used to analyze various shell problems. 
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INCREMENTAL. TQTAL-LAGRANGIAN FORMULATION OF A CONTINUOUS MEDIUM 

The primary objective of this section Is to review the formulation 
of equations governing geometrically nonlinear motion of a continuous 
medium. In the Interest of brevity }nly necessary equations are 
presented. For additional details the reader Is referred to References 
[ 6 - 10 ]. 

We describe the motion of a continuous body In a cartesian 
coordinate system. The sirultanecus position of all material points 
(l.e., the configuration) of the body at time t Is denoted by C t , 
and C Q and denote the configurations at reference time t = 0 and 

time t + at, respectively (see Fig. 1). In the updated Lagranglan 
description all kinetic and kinematic variables are referred to the 
current configuration at each time and load step. In the total 
Lagranglan description all dependent variables are referred to the 
reference configuration. The updated Lagranglan is more suitable for 
motions that Involve very large distortions of the body (e.g., high- 
velocity Impact). The total Lagranglan Is more convenient for motions 
that Involve only moderately large deformations. In the present study 
the total Lagranglan formulation Is adopted. 

Here we present a derivation of the equilibrium equations at 
different time steps using the total Lagranglan approach. The 
coordinates of a typical point In C t Is denoted by t x = ( fc Xj , t x ? , t x 3 ) . 
The displacement of a particle at time t Is given by 


\ = t x - °x or t u j = t x 1 -°x 1 (1) 

The Increment of displacement during time t to t + At Is defined by 


t+at 


U 1 - U 1 


( 2 ) 


I 



Figure 
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The principle of virtual displacements can be employed to write tie 
equilibrium equations at any fixed time t. The principle, applied to 


the large-displacements case, can be expressed mathematically as 

t+At; i A\t * f t+At< 


; 


•j, dv o + ; 


" S 1J s < t+ %> dV o 


J 


t+At 


T. 5U. dA„ + r 

1 1 0 


»o tM ‘ F 1 SU 1 dV o 


(3) 


where summation on repeated indices is implied; V 0 , A Q , and p Q denote, 
respectively, a volume element, area element, and density In the Initial 
configuration, S^j are the components of second Piola-KIrchhoff stress 
tensor, e,, the components of Green-Lagrangian strai n tensor , the 

components of boundary stresses, and are the components of the body 
force vector; the superi • >ed dots on u^ denotes differentiation with 
respect to t^me, and 6 Jenotes the variational symbol. In writing Eq. 
(3) it Is assumed that is related to the displacement components by 
the kinematic relations 


t+At 




1 

2 


( t + At u 


1J 


t+At, 


t+At, 


t+At 


j,1 


m.i 


Vj> 


(4) 


*. L A *• 

where j = au^/axj. The strain components e jj can be expressed in 
terms of current strain and incremental strain components as 


t+At _ 1 ,t , t , t.. t, , 

1 J 2 1 , j j , i m,i 

1 t t 1 

+ 2 ^ U i,j + U j,i + U m,1 U !Ti,j + u m,i u m,j^ + 2 u m,1 u m,j 

i * (e,j + n,j) 


(5) 
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where e^j and denote the linear and nonlinear incremental strains. 

4. 11 |i 

The stress components S^j can be decomposed Into two parts: 

1 * ts ij = t$ 1j + S 1j 

where S^j is the incremental stress tensor. The incremental stress 
components S-jj are related to the incremental Green-Leg^ange strain 
components, = e^ + by the generalized Hooke's law: 

S ij = C, i jka e k£* 

where are the components of the elasticity tensor. Using Eq. (4)- 

(7), Eq. (3) can be expressed In the alternate form 


L °o t+4t u 1 4u 1 dV o + L C tjki (e kn 5,| 1j * "k« 5e 1j> dV o 
o u o 

+ / ‘s,j dV 0 - <M - / ‘s,j «„ (8) 

o o 

where sW Is the virtual work due to external loads, 

FINITE-ELEMENT MODEL 

Geometry of the E lement 

Consider the solid three-dimensional element shown in Fig, 2. 

The coordinates of a typical point In the element can be written as 


x., = 


n 

L 

j=l 


*j< 5 


ja? 2 ) 2 


1±£ 


(x}) 


top 






) 1^. ( X J\ 

' 2 bottom 


(9) 


where n Is the number of nodes, are the finite-element 

interpolation (or shape) functions, which take in the element, the value 
of unity at node 1 and zero at all other nodes, ^ and are the 
normalized curvilinear coordinates in the middle plane of the shell, 
and c is a linear coordinate in the thickness direction and xj, x^, and 
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are the global coordinates at node 1. Here ?j S c 2 , and z are assumed 
to vary between -1 and +1. Now let (see Fig. 2) 



^Vtop 


- K) 


i, 

k bottom 


V , v 1 / 

-3 * - 3 /l 




( 10 ) 


wnere v^ Is the k-th component of the vector v^. Then Eg. (9) becomes 

X 1 * j, + *jf h j *3i J C'D 

where hj Is the thickness of the element at node j. For small 
deformation, the displacement of every point In the element can be 
written as 


U i = jf 1 *j Iu 1 + 5 I Hl 9 2 ‘ e 21 e l )] 02) 

1 i * i *• 1 

where e.j and aJ, are the rotations about (local) unit vectors e^ and e 2 » 

respectively, u-j, u 2 , and U 3 are the displacement components 

corresponding to the global coordinates x l9 x 2 , x 3 directions 

respectively, and uj 9 u^ and u^ are the values of the displacements 

(referred to x) at node 1. In writing Eq. (12), we assumed that a line 

that is straight and normal to the middle surface before deformation is 

still straight but not necessarily 'normal' to the middle surface after 

deformation. The strain energy corresponding to stress perpendicular to 

the middle surface Is Ignored to Improve numerical conditioning when the 

three dimensional element is employed. This constraint corresponds only 
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to a part of the usual assumptions of a two-dimensional shell theory. 
The relaxation of the requirement that straight lines perpendicular to 
the middle surface remain normal to the deformed middle surface permits 
the shell to experience shear deformation - an important feature in 
thick shell situations. 

Displacement Field In the Element 

In the present study the current coordinates ^ are Interpolated 
by the expression 


% • f, ’/‘I + 7 ‘"SI 1 


(13) 


and the displacement by 

* U 1 • jf, *J ftg i + i' h j (* e 3l - < 14 > 

U 1 * J f 1 *jl“l + \ ( t+4t =31 - (15) 

Here ^u^ and denote, respectively, the displacement and incremental 
displacement components In the x^-dlrectlon at the j-th node. The unit 
vectors ej and e,j, can be obtained from the relations 


e{ = (E 2 x t e^)/!E 2 x 


*1 t - ■? tM 

?2 = h* Si 


(16) 


where is the unit vector along the (global) X 2 ~axis. If we assume 
that the angles ej and are very small, then we can write 
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?3 


- ‘iM ♦ ‘i!4 


( 17 ) 


Substituting Eq. (17) Into Eq. (15), we obtain 


n 

z 

j=l 


*j Iu 1 


1 

2 


Ch j ( ' 


t:j 9 J 

e 2i 9 1 




e ii s i)i 


(18a) 


or 

{u} = IT] {a} (18b) 

where {u} Is the column of three displacements at a point, {a} Is the 
column of 5n (five per node) displacements: u^, Q'j, aij, j = l,2,...,n; 1 
= 1,2,3, and [T] Is the transformation matrix defined by Eq. (18a). 

Thus for each time step one can find the normal vectors from Eq. (16) 
and (17), and the incremental displacements at each point from Eq. (18) 
once the five generalized displacements at each node are known. 


Element Stiffness Matrix 


The strain-displacement equations (4) can be expressed in the 
operator form 

{e} - (A]{u 0 } (19) 

where {e} = {e^ e 22 e 33 2e^ 2 2e 13 2e 23 J T , [A] is a function 

of *u . and (u } is the vector of the components of the displacement 
o I , j o 

gradient 

(u Q } = (\i U- | s 2 U l,3 U 2,l U 2,2 U 2,3 U 3,l U 3,2 U 3 S 3^ (20) 

The vectors {u } and {e} are related to the displacement increments by 


{u Q } = [N]{u} = [N] [T] {a} 
{e} = [A] [ N J [T] {a} = [B]{a} 


( 21 ) 

( 22 ) 


where [ N ] is the operator of differentials. 


Substitution of Eq. (22) into Eq. (8) yields 


; 0 o [T] t (i}dV o + (‘[K l J + t lK Nl l)(a) - t+at {fi) - t+4t {Fl (?3) 
v o 

t* f* 

where [KJ, (K N) J» {R}, and {F} are the linear and nonlinear stiffness 
matrices, force vector, and unbalanced force vectors; 

t [K L ] = f t [B] T [C] 1: fBjdV 0 , t (K NL ) - f t f 8 ] T ( S ] t [S]dV 0 
V o V o 

{F} = J t lB) T {S}dV (24) 

V o 

Here [Sj and {$} denote the matrix and vector, respectively, of the 
second Piola-KIrchhoff stress. 

Since we are dealing with laminated composite structures, the 
important thing Is how to perform the Integration through the 
thickness. One way Is to pick Gaussian points through the thickness 
direction. This Increases the computational time as the number of 
layers Is Increased, because the integration should be performed 
separately for each layer. An alternative way Is to perform explicit 
integration through the thickness and reduce the problem to a two 
dimensional one. The Jacobian matrix, In general, Is a function 
of 4nd t* The terms In $ to the first power may be neglected, 

provided the thickness to curvature ratios are small. This 
approximation implies that derivative of x^ with respect 
to Sj, and c are substantially the same at either end of a mid- 
surface-normal line. Thus the Jacobian (J) becomes independent of e and 
explicit integration can be employed. If 5 terms are retained In [JJ, 
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i 

I 

| 

Gaussian points through the thickness should be added. In the present i 

study. It Is assumed that the Jacobian Is Independent of 5 . ] 

* | 

Time Integration ) 

The Newmark Integration scheme Is used to convert the ordinary ; 

differential equations in time, Eq. (23), to algebraic equations. In l 

t- 

the Newmark scheme, displacements and accelerations are approximated by ' 

t 

t+lt {a} = ‘{4} + 4t 2 (i} + [(i - 8)‘{i} ♦ B t+it {a}I<at) r { 

l 

i 

» ‘(i) + [(1 - T )‘{i} + r t+at (i}l4t (25) 1 

•i 

where {a} Is the generalized displacement vector of any point 
and a and Y are the dimensionless parameters of the approximation. For 
the constant average acceleration case, we have s = and y = and for 
the linear acceleration method a = ^ and Y = | (see [ 11 ]). 

Substituting Eq. (25) Into Eq. (23), and some algebraic 
manipulation leads to 

(a Q t [MJ + ‘[K]){ 4 < k >} - t+ 4 t (R} - t+At {F (k-1) } + 

+ a 2 (‘{Pi) - -‘( p 3 ))l (26) 

where 

a ° 3 Twi? 9 a? = ^ 9 a 3 = - 1 * and 

[MI = J P n ^[T] T t [T] f ; 
v 0 


i 
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( p l> ■ j„ »0 ‘W t T l dv c 


fp 2 } - ; ° 0 l+st (4} (k ' 1) [T]dv ( 

V o 

t P 3 l ° * f » 0 ‘WITIdV 0 


( p al ■ ! »„ ‘{i} (TldV 


4 J J v "o 
0 


(27) 


This completes the finite-element formulation of the 3-D degenerated 
element. 


DISCUSSION OF THE NUMERICAL RESULTS 

The results to be discussed are grouped Into three major 
categories: (1) static bending, (2) natural vibration, and (3) 
transient response. All results, except for the vibrations, are 
presented In a graphical form. All of the results presented here were 
obtained on an IBM 370/3081 computer with double precision arlthmatlc. 

Static Analysis 

Here we present a discussion of four example problems, all 
Involving shell structures. 

1. Cylindrical Shell Subjected to Radial Pressure Consider a 
circular cyllnd-ical panel of the type shown In Fig. 3. The shell is 
clamped along all four edges and subjected to uniform radial Inward 
pressure. The loading is nonconservative, that is, the direction of the 
applied load Is normal to the cylindrical surface at any time during the 
deformation. The geometric and material properties are 


R 


free edge 


supported 
by a rigid 

diaphram: u 


Figure 3 



Geometry of the cylindrical shell used in Problem 1 
of the static analysis. 
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R = 2540 mm, a = b = 254 (ran, h a 3.175 mm, 

0 = 0.1 rad, £ 13 3.10275 kN/mm 2 , v = 0.3 
Due to the symmetry of the geometry and deformation, only one quarter of 
the panel Is analyzed. A load step of 0.5 KN/m^ was used In order to 
get a close representation of the deformation path. Fig. 4 contains the 
plot of central deflection versus the pressure. The solution agrees 
very closely with that obtained by Ohatt [13]. 

2. Orthotropic Cylinder Subjected to Internal Pressure Consider a 
clamped orthotropic (£2 = 20 x 10 6 psl, £^2 = 3.75, G 12 /E2 = 

0.625, v = 0.25) cylinder of radius R = 20" and length 20", and 
subjected to Internal pressure. p Q = 6,41 /w psl, A mesh of 2x2 nine- 
node elements Is used to analyze the problem. The linear center 
deflections obtalneo by the 2-D and 3-D elements are 0.0003764 In., and 
0.0003739 In., respectively. These values compare favorably with 
0.000366 In. of Rao [14] and 0.000367 in. of Timoshenko's analytical 
solution [15]. The latter two solutions are based on the classical 
shell theory. 

In the large-deflection analysis the present results are compared 
with those of Reference 5. A value of 2.5 ksl is used for the load 
step. Figure 5 contains a comparison of the present deflection with 
that of Reference 5, which used a 3-D degenerated element based on the 
updated Lagranglan approach. The agreement is very good. 

3. Nine-Layer Cross-Ply (oygoVOV. . .) Spherical Shell Subjected 

to Uniform Loading Consider a spherical shell laminated of nine layers 
of graphite-epoxy material ( E -| /H 2 = 40, ^ 12/^2 = ^13 = G 12 “ 

g 23* v I 2 ~' 25 sub J ected t0 uniformly distributed loading, and simply 
supported on all its edges (i.e., transverse deflection and tangential 


\ 




6Q 



Figure 5 Center transverse deflection versus 
internal pressure 
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rotations are zero'. A comparison of the load-deflection curves 
obtained by the present elements with those obtained by Noor [4J is 
presented (for the parameters h/a => 0.01 and R/a = 10) in Fig. 6. The 
results agree very well with each other, the present 2-D results being 
closer to Noor's solution. This is expected because Noor's element is 
based on a shell theory. 

4 . Two-Layer Cross-Ply a n d Angle-Ply (45°/-45°) Shells Under 
U niform loading The geometry of the cylindrical shell used here is the 
same as that shown In Fig. 3. The shell is assumed to be simply 
supported on all edges. The material properties of individual lamina 
are the same as those used in Problem 3. A mesh of 2x2 nine-node 
elements in a quarter shell is used to model the problem. The results 
of the analysis are presented in the form of load -deflection curves in 
Fig. 7. From the results, one can conclude that the angle-ply shell i ,: . 
more stiffer than the cross-ply shell. 

The geometry and boundary conditions used for the spherical she 1 Is 
are the same as those used in Problem 3. The geometric parameters used 
are: R/a - 10, a/h = 100. The load-deflection curves for the cross-ply 

and angle-, 1,y shells are shown in Fig. 8. From the plot it is appa-ent 
that, for the load range considered, the angle-ply shell, being stiffer, 
does not exhibit much geometric nonlinearity. The load-deflection curve 
of the cross-ply shell exhibits varying degree of nonlinearity with the 
load. For load values between 100 and 150, the shell becomes relatively 
more flexible. 

Natural Vibration of Cantilevered Twisted Plates 

Here we discuss the results obtained for natural frequencies of 
various twisted plates. This analysis was motivated by their relevance 


Center deflection, w/h 



Figure 6 Deflection versus load parameter for a nine- 
layere cross-ply (0°/90°/0 0 / . . . ) spherical 
shell 



HondiPiensional ized center deflection, w/h 
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to natural vibrations of turbine blades. Consider an Isotropic 
cylindrical panel with a twist angle 0 at the free end. Table 1 
contains the natural frequencies of a square plate for various values of 
the twist angle e and ratios of side to thickness. A 2x2 mesh and 4x4 
mesh of 9-node elements are employed to study the convergence trend. 

The results of the refined mesh are Included In the parentheses. The 
results obtained by using the 4x4 mesh are lower than those predicted by 
the 2x2 mesh, showing the convergence. The results agree with many 
others published in a recent NASA report. Table 2 contains natural 
frequencies of twisted plates for the aspect ratio of 3. 

Transient Analysis 

1. Spherical Cap Under Axisymmetric Pressure Loading Consider a 
spherical cap, clamped on the boundary and subjected to axisymmetric 
pressure loading, p Q . The geometric and material properties are 

R = 22.27 in., h = 0.41 in., E = 10.5 x 10® psi, \> = 0.3, 
p = 0.095 lb/1n 3 , 0 = 26.67°, p Q = 100 psi. At = 10~ 5 sec. 

This problem has been analyzed by Stricklin, et al. [16] using an 
axlsymmettic shell element. In the present study the spherical cap is 
discretized into five nine-node 2-0 and 3-0 elements. Figure 9 contains 
the plot of center deflection versus time. The present solutions 
obtained using the 3-0 and 2-0 elements are in excellent agreement in 
most places with that of Stricklin et al [16]. The difference between 
the solutions is mostly in the regions of local minimum and maximum. 

2. Two-Layer Cross-Ply Plate Under Uniform Load A cylindrical 

shell with a = b = 5", R = 10 1 ', h = 0.1" is simply-supported on the four 

edges. Is analyzed. The shell Is laminated by 2 layers (0°/90°) and 

- a 4 P 

exerted by a uniform step load P = j = 50. r gure 10 contains a plot 
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Table 


Natural Frequencies of Twi sted Square Plates 


aa Z /oR7D , 0 = , v = 0.3 

1 2(1 -v ) 


a 

h 

Twist 
Angl e 

1 

2 

Mode 

3 4 

5 

6 


0° 

* 3.4556 

+ (3.4583) 

8.41 10 
(3.3353) 

22.0999 

(21.0238) 

28.2089 
(26. 7465) 

31 .9740 
(30.1454) 

55.1625 

(52.0704', 


15° 

3.4359 

10.2920 

21 .5199 

27.2054 

32.7430 

44.5375 

20 

30- 

3.3790 

(3.3694) 

13.7014 

(14.2222) 

19.9840 

(18.9795) 

25,0943 

(26.8104) 

34.3341 

(34.4591) 

45.8987 

(45.7547) 


45° 

3.2908 

13.1009 

15.9097 

23.5680 

35.5332 

45.7013 


60° 

6.1800 

17.8319 

15.5635 

24.1842 

36.1466 

44.9152 


0° 

* 3 .33916 

**(3.3390) 

7.3948 
(7 .3559) 

10.8083 

(10.883) 

13.4930 
(1 7.757) 

23.7907 

(22.769) 

26 .0552 
(24.125) 


1 5° 

3.31713 

(3.3170) 

7.4816 

(7.4504 

10.8053 

(10.774) 

18.4043 

(17.771) 

23.6767 

(22.694) 

24.9474 

(24.083) 

5 

30° 

3 .2538 
(3.2538) 

7.7593 

(7.7089) 

10.5248 

(10.478) 

18.4091 

(17.795) 

23 .3734 
(22.471) 

24.61 16 
(23.943) 


45° 

3.1570 
(3.1 569) 

8.1435 

(8.0728) 

10.1270 

(10.062) 

18.3843 

(17.79) 

22.9126 
(22.1 17) 

24.0566 
(23.651 ) 


60° 

3 .0370 
(3. 0366) 

8.5855 

(8.4814) 

9.67198 

(8.5911) 

18.3089 
(1 7.730) 

22.3670 

(21.684) 

23 .3533 
(23.160) 


* 2x2, 9~node mesh 
**3x3, 9 - node mesh 
+ 4x4, 9-node mesh 
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Table 2 Natural Frequencies of Twisted Rectangular Plates 
(b/a = 3, 3x3 mesh of nine-node elements) 




w = <jjb 

r/’oTTU n = E 

IF . = 

0.3 



12(1 

? ' 
!-v) 2 

a 

Tvyi s t 



Made 




h 

Angle 

1 

2 3 

4 

5 

6 

7 


0° 

3.4150 

20.8772 21.6190 

65.9706 

66.2590 


127 .256 


15° 

3.4009 

20.8798 22.1118 

21.5032 

68.0938 

69.3253 

130.284 

20 

30° 

1 .3598 

19.4048 25.3743 

60.2183 

73.5180 

77 .4493 

138.176 


45° 

3.2956 

17.5289 29.8404 

58.2600 

80.9488 

88.5246 

148.8975 


60° 

3.2136 

1 5.7431 34.8827 

55.8921 

89.2028 

100.7760 

155.070 


6° 

3.3908 

15.551 

19.124 

21.065 

59.924 

61.949 


15° 

3.3161 

’8.192 

19.231 

21 .572 

60.088 

60.830 

5 

30° 

3.3336 

14.379 

19.549 

22.81 1 

60.576 

58.472 


45° 

3 .2674 

13.449 

20.060 

24.404 

61 .360 

55.374 


60° 

3.1833 

12.548 

20. 741 

26.139 

62.416 

53.381 



Figure 9 


Center transverse displacement versus time 
for a spherical cap under axi symmetric dynami 
loading ( load = 100 psi.) 





Center deflection, wxlO 


3-0 Element 
2-0 Element 


Time, txlO 4 

Figure 10 Center deflection versus time for t 
cress-ply cylindrical shell subject 
uniform step load 
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of the center deflection versus time for 2-0 and 3-D elements. The time 
step used Is at = 0.1 x 10 _<: * sec. The solutions obtained using the two 
elements are In good agreement. 

3. Two-Layer Angle-Ply (45°/-45°) Spherical Shell Under Uniform 
Loading Consider a spherical shell with a = b = 10" „ R = 20" and h = 
0.1", simply supported at four edges and Is exerted by a uniform step 
load. The shell consists of two layers, (45V-45 9 ). Figure 11 contains 
the plot of center deflection versus time for P = 50 and P = 500 with 
time step 0.2 x 10”^ sec. For the small load the curve Is relatively 
smooth compared to that of the larger load. This is due to the fact 
that the geometric nonlinearity exhibited at P = 50 Is smaller compared 
to that at ? = 500. 

CONCLUSIONS 

The present 3-0 degenerated element has computational simplicity 
over a fully three-dimensional element, such as those developed in [17], 
and the element accounts for full geometric nonlinearities in contrast 
to 2-D elements based on shell theories. As demonstrated via numerical 
examples, the deflections obtained by the 2-D shell element deviate from 
those obtained by the 3-D element for deep shells. Further, the 3-D 
element can be used to model general shells that are not necessarily 
doubly-curved. For example, the vibration of twisted plates cannot be 
studied using the 2-D shell element discussed in [12]. Of course, the 
3-D degenerated element is computationally more demanding than the 2-0 
shell theory element for a given problem. In summary, the present 3-D 
element is an efficient element for the analysis of laminated composite 
plates and shells undergoing large displacements and transient motion. 

The 3-D element presented herein can be modified to Include thermal 
stress analysis capability and material non! Inearl ti While the 
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Inclusion of thermal stresses is a simple exercise, the inclusion of 
nonlinear material effects Is a difficult task (see (18-201). An 
acceptable material model should be a general izatfon of Ramberg-Osgood 
relation to an anisotropic medium. Another area that requires further 
study is the inclusion of damping effects, which are more significant 
than the shear deformation effects. 
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PART 3 


NONLINEAR MATERIAL MODELS FOR COMPOSITE PLATES AND SHELLS 

K. Chandrashekhara and J. N. Reddy 
Department of Engineering Science and Mechanics 


SUMMARY 

Nonlinear material models for laminated structures are described 
and their Incorporation In the finite-element formulation of laminated 
plates and shells Is presented. Numerical results for several sample 
problems plates and shells are presented and validated by comparison 
with those available In the literature. 


INTRO DUCTION 

Composite materials are known to exhibit significant non- 
linearities In stress-strain behaviour even at low strains. Most of the 
currently used matrix materials In composites have high strain 
capabilities and the investigation of the bending of composite shells 
undergoing large deformation, yielding is apt to occur and its effect 
mu^t be accounted for in the analysis. The nonlinearity is not 
isotropic but varies with direction, as do the elastic properties. 

Models for such elastic-plastic behavior of orthotropic and anisotropic 
materials are not well developed. 

The total stress-strain laws are mathematically more convenient 
than Incremental laws but are physically not sound. The criterion 
approximately describing the yielding of Isotropic material is that of 
von-Mises. The simplest yield criterion for anisotropic material is 
therefore one which reduces to von-Mises law when the anisotropy is 
vanishingly small. Hill's yield criteria assumes relatively sli.ple ease 
of orrhotroplc anisotropy, that Is, there are three mutually orthogonal 
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planes of symmetry at every point and the intersection of these planes 
are concldered as the principal axes of anisotropy. Fiber reinforced 
composite structures almost invariably possess this kind of symmetry. 

In the present study a nonlinear material model Is developed for 
composite plates and shells, and numerical results for bending are 
presented using the finite element method as exact solutions are not 
tractable for elastic-plastic p. blems involving complex geometries. 

MATERIAL MOD EL 

In the present model. Hill's anisotropic yield criteria for 
elsatlc-perfectly- plastic material is used. Hill's (1] yield function 
Is, 

f(o^j) = F(o 2 - a 3 ) 2 + G ( a 3 - a ^) 2 + H(oj - a ,,) 2 

+ 2La 2 3 + 2 Mo 2 3 + 2Na 2 2 =1 (1) 

where F, G, H, L, M, N are parameters characteristic of the current 
state of anisotropy given by. 



and X, Y, Z are the tensile yield stresses in the principal direction of 
anisotropy and R, S, T, are the yield stresses in shear with respect to 
the principal axes of anisotropy. 

It should be noted that Hill's criteria is based on the assumption 
that the superposition of a hydrostatic stress does not influence 
yielding and there is no Bauschinger effect. Also, the yield criterion 
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has this form when the principal axes of anisotropy are the axes of 
references. 

For a plane stress state in the 1-2 plane with transverse shear,, 
equation (1) reduces to: 

f = (G + H)o* + (F + H)ajj - 2Ha^2 

+ 21 ofj 3 + 2Ma\ 3 + 2Hc\ 2 = 1 (2) 

For an isotropic material:. 


X = Y = Z = o Q , 

the yield stress in uniaxial ' 'nslon and according to the vDn-MIses 
yield criteria 12] 


R = $ = T = 

/3 

Therefore, F = C = H = and 2L = 2M = 2N = ^ and equation (2) 

2a a 

becomes , 0 

f . 2 A 2 > 2 A 2 ^ 2 , 2 

f - °i + °2 ‘ a l CT 2 + 3(a 23 + a 13 + G W 3 u o 

which Is the familiar von-MIses yield criteria- 

If the principal axes of anisotropy 1,2 do not coincide withthe 

reference axes x, y, but are rotated by an angle a, then the stresses In 

equation (2) are obtained using the transformation as: 

2 2 

a-i^a cos g+o sin a + a sina cose 
lx y xy 

2 2 

o- « a sin a + a cos a - a sine cos 8 

c a y Ay 

°p? = " a v 7 sino + a cosa 
°i3 “ °xz “ s9 + V s1ne 

2 ? 

o, 0 = -2a sine cose + 2cr sine cose + a vl (cos e - sin a) 
x y xy 


Elastic-Plastic Constitutive Equations 

In the incremental theory of plasticity, the total strain Increment 
is the sum of the elastic and plastic components 
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de = de e + de P (3) 

The elastic strain Increment is related to the stress increment by 
Hooke's law as, 

dr e = [D e ]~*dcr (4) 

where [D e J is the elastic modulus matrix which for orthotropic material 
takes the form, 


[B 


ei = 


1 


1 ”’ u 12 v 21 

V 12 E 2 

1 -' , 12 v 21 

0 

0 

0 


v 12 E 2 
1 ■" v , 


0 


12 21 
En 


1 ’ V 1Z V 21 


0 

0 

0 


a 23 

0 


0 

0 

0 

d 13 

0 


0 

0 

0 

0 

a 12 


( 5 ) 


The normality rule for an associated plastic flm. is. 


de p = d> 


af 

3a 


where dx Is the positive proportional Ity constant, evaluated using the 
condition that during the plastic deformation, the stresses remain on 
the yield surface so that. 


df = If do = 0 
3o 


The stress-strain relation in the plastic range is given by 1 3 ] 

da = [D ep ]de 

where 


[D ep ] - [0 e ] - - 


!“ e H£}^) T ID e l 


{— } T (D e ] {— 1 
l 3o ; 1 u 3a ; 


( 6 ) 
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Hence the modification called for In the elastic-plastic analysis 

would be solely the replacement of the elasticity matrix [D e ] by the 

elastic-plastic matrix (D ep ] for the yielded elements at the successive 

stages of calculation. It should be noted that the [C ep ] matrix Is 

populated and accordingly the transformation of the stress-strain 

relation from the material axes, {{a} = to the shell 

coordinate axes, ({a} = [D ep l{e}) 9 will be modified as shown In 

xy 

Appendix I. 

FINITE ELEMENT FORMULATIO N 

Consider a laminated shell constructed of a finite number of 
uniform thickness orthctropic layers, oriented arbitrarily with respect 
to the shell coordinates Uj,5 2 ,c). The orthogonal curvilinear 
coordinate system Is chosen such that and curves are 

lines of curvature on the midsurface ;=0, and ^-curves are straight 
lines perpendicular to the surface r,=0. 

For the small displacement Sanders shell theory which accounts for 
transverse shear deformation, the strain displacement relations are 
given by { 5 ] , 

e 1 = e i + CK 1 

o . !*1 

e l ax x Rj_ 9 K 1 " axj 

o _ au 2 U 3 . 3 *2 

e 2 3X 1 R 2 9 k 2 ~ ax 2 

_ 3U. 3U 0 3tt» 1 3 3U- 3U, 

- ~ . x _ p / 

'6 3X 2 3Xj 9 6 3X 2 3Xj O V 3X^ 3X 2 7 

o 4 «Z 

e 4 ■ *2 3X 2 ' R 2 


where 


30 


o 

e 5 3 


3U. 
’1 + 3^ 


1 ,L 

2 { R. 


d •. 


« 1 de 1 


(i - 1,2) 


Here (1 = 1,2) are the principal radii of curvature, u^ are the 
displacements of the reference surface along = 0 axes, and <t> 2 

are the rotations of the transverse normals about the and Cj-axes 
respectively. 

The stress-strain relations, transformed to the shell coordinates, 
are of the form 

M = [Q] {e } 

where - are the material properties of k th -layer (see Append-* ■ I). 

The principle of virtual work for the present problem Is given by 


0 * E 5e g + 


+ «e g * q6U 3 }o 1 a 2 d^^dC 2 ]dc 


(7a) 


= j" [N^Se^ + + N g Se g + + MgOKg 

+ ^4 5e 4 + ^5 fie 5 " qSUjJo^o^dS jjd£ 2 (^b) 

where q is the distributed transverse load, and are the stresses 
and moment resultants. 


L 'k 

(» i *M 1 ) = I J o^l.cjdc (1 = 1,2, 6, 4, 5) 
k=:1 c k-l 


Here (c^ ^, 5 ^} are the t-coordi nates of the k th layer, and L Is the 
total number of layers In the laminated shell. 
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It should be noted that the equations of equilibrium can be derived 
from Eq. (7b) by Integrating the displacement gradients In e° by parts 
and setting the coefficients of su 1 to zero separately. We obtain 

1 + -i- (N fi + c„M e ) + ^ - 0 


ax 


1 


ax. 


'o”6 J 

aN, 


1 


N, 


TxJ’ (N 6 " c o M 6 ) + a^ + Rj = 0 
aN,- aN. N, N ? 
3x^ + ax^' ( R7 + Rj" q) = ° 
aM, aM fi 

1x7 + IxZ " N 5 = 0 

aM, aM 

^ + a4" N4 = ° 

The resultants (N.j, M^) are related to (e^,^) by, 
o 

+ i5„„tc„ 

1.J = 1,2,6, 4, 5 


N 4 = ^ij E j + 8„_>c 


'1 


up P 


= 8 .e'l + D * 

& aj J ap P 


a,p = 1 ,2 V 6 (with i=1 for 1 = 1,2,6) 


(8) 


Here A^j, 3^ and D^j denote the oxtenslonal , flexura'i-extenslonal 
coupling, and flexural stiffnesses of the laminate: 


<A ij* B ij*°ij ) ' A ! Q i ( j > 

K-i 5 j< _ 1 


( 9 ) 


In the unabridged notation equation (8) takes the form: 
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' N 1 


A 11 

A 12 

A 16 

-14 

^15 

d ll 

B 12 

i 

U> 

DO 


o a 
e l 

n 2 


A 12 

A 22 

A 26 

-24 

A 25 

B 12 

B ?2 

B 26 


0 

e 2 

n 6 


A 16 

A 26 

A 66 

^46 

-56 

B 16 

B 26 

b 66 


0 

e 6 

n 4 


-14 

-24 

-46 

A 44 

A 45 

-14 

-24 

-46 


0 

e 4 


‘3 











M S 


-15 

-25 

-56 

A 45 

A 55 

-15 

-25 

-56 

a 

0 

e 5 

Ml 


8 U 

8 12 

8 16 

-14 

^15 

D 11 

°12 

°16 


r l 

m 2 


B 12 

B 22 

B 26 

-24 

-25 

D 11 

D ?2 

°26 


*2 

M 6. 



B 26 

B 66 

-46 

-56 

°16 

°2Q 

°66_ 


*6 

^ / 


The underscored coefficients are due to material nonlinear stress-strain 
relationship. It should be noted that the coefficients A 44 , A 45 and Ag 5 
defined 1r iquation (3) has to be corrected for the parabolic variation 
of the transverse shear stress, as 

L ? k 

^ A 44 sA 45® A 55) = E J* ^1^44^* k l k 2^45^ k 2^55^ dc 

k=1 Vl 

where k are the shear correction factor. 

A typical finite element Is a doubly-curved shell element whose 
projection Is an Isoparametric rectangular element. Over the typical 
shell the displacements a^e Interpolated by 

expressions of the form. 


N 


U 1 = / u 1^j^ x l ,x 2^ 8 1 = 

j — 

N 1 

* 1 35 £ *4<> 4 (x 1( x 9 ) , 1-1,2 


( 12 ) 


j=l 


vj v 1* 2 


where tpj are the interpolation functions, and u:j and are the nodal 
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values of u^ and 4 > 1S respectively. For a nine node quadratic element 
the element stiffness matrix Is of order 45x45. 

Substitution of equation (12) Into the virtual work principle, Eq. 
(7b) yields an element equation of the form 

[Kl{4} = {F} (13) 

where { 4 } = {{uj}, {u 2 J, {u 3 J, {♦j}, {$ 2 }} T , (K] Is the element 
stiffness and {F} Is the forco sector. In t^e Interest of brevity, the 
coefficients of stiffness matrices are Included 1r; Appendix II. 

It should ue noted that the underscored coefficients In Eq. (10) 
are also redefined like the shear coefficients In Eq. (11) and reduced 
Integration is performed for the terms arising In the element stiffness 
matrices due to the presence of these coeffi:1ents to avoid the so- 
called locking effect. 


NUMERICAL RESULTS 
The Parameters of Anisotropy 

When considering the modeling of a material system, one must always 
survey the availability of material property data. In the present 
theory, to describe fully the state of anisotropy, the six Independent 
yield stresses in Hill's criteria are needed to be ki’own from uniaxial 
tests. For numerical results, two typical composite materials namely, 
boron/epoxy and graphite/epoxy are considered with the following 
material constants: 

Boron/Epoxy 

Ej = 30.0 x 10 6 psl , E 2 = 3.2 x 10 6 psi 
G^ 2 « 1.05 x 10® psl , v^ 2 = 0.21 , G 23 = Gj 3 = G^ 2 

X = 195 x 10 3 ps1 ; Y = l = 12.5 x 10 3 psl 
R = S => T = 18.0 x 10 3 psi 
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Graphite/Epoxy 

Ej - 18.88 x 10 6 psl ; E ? = 1.376 x 10 6 psi 
Gjg 3 0.688 x 10® psl ; 3 0.343 ; = Gjj a 0^ 

X = 222.7 x 10 3 psl ; Y •-» Z a 6.35 x 10 3 psl 
R =- S = T - 9.92 x 10 3 psi 


Solution Procedure 

The solution of the elastic plastic problem Is reached by an 
incremental and iterative procedure. The direct iteration technique is 
followed in the present analysis. 

For each load increment, the system of equations are established by 
assembly! r.g the element matrices and the displacement {a} Is obtained 
from Eq. (13) . Consequently, the state of stress and the value 
of f (cr^ j) are calculated for each element If f < 0, then the process 
Is elastic and the material matrix Is obtained from equation (5). If f 
> 0, then the total stresses are readjusted so as to make f a 0 and the 
elastic-plastic matrix is calculated from Eq. (6). Once the convergence 
is achieved, the next load Increment is applied and the iteration 
procedure is repeated. 

If the application of a small load Increment causes very large 
deflection, the calculation is stopped and the limit lead is considered 
to be found. 

Sample Problems 

The present elastic-perfect ly plastic formulation Is applied to a 
variety of bending problems using 2x2 mesh of a nine noded quadratic 
element. The shear correction factors k^ = k^ were taken to be 5/6. 
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All computations were made using an ICM 3081 processor with double 
precision arithmetic. 

The results of the sample problems are presented and compared. If 
possible, with the existing solutions to evaluate the present 
formulation. 

1. Cylindrical Shell Roof A cylindrical shell subjected to 
uniform vertical loading Is considered. Due to symmetry, only a 
quadrant of the shell was analyzed. The geometry and modeling of the 
shell roof are shown In Fig. 1. The material behaviour Is studied with 
the properties: 

E l = E 2 = 2.1 x 10 4 MN/m 2 ; v = 0.0; 

G 12 = 1.05 x 10 4 MN/m 2 ; G 23 = G 13 = G 12 
X = Y = Z = 4.1 Mn/m 2 ; R 3 S :i T = 2.367 MN/m 2 

The results obtained for the vertical displacement at the central 
point of the free edge A versus loading was shown in Fig. 1. The 
solution obtained compares well with those reported In Ref. [6]. The 
apparent discrepancy can possibly be due to a different boundary 
condition on the curved edges and the type of material model used. 

2. Simply-Supported Square Plate A uniformly loaded simply 
supported square plate was studied in the second example. The geometry 
of the plate Is shown in Fig. 2. The following material properties were 
considered: 

E x = E 2 = 10 x 10 6 psi ; v = 0.3 
G 12 = 3.846 x 10 6 psi ; G 23 = G 13 = G 12 
X = Y = Z = 144,000 psi ; R - S = T = 83,138.4 psi 





Figure 2. 


Present solution 



| 1 1 r 

.02 .03 .04 .05 

9 

Center deflection, wD/(4M o a ) 


Load-deflection curves for a simply-supported 
square plate under uniform transverse load 
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A non-dimensional ized plot of the centre displacement of the plate 
versus the load are shown In Fig. 2. The results are compared with 
those presented In Ref. [7]. 

3. Two Layer Cross-Ply [0/90] and Angle-Ply (-45/45) Simply 
Supported Spherical Shells Figure 3 contains the results for the cross- 
ply shell made of two typical materials, namely, boron/epoxy and 
graphite/epoxy under uniform load. For a given load, the shell made of 
graphite/epoxy deflects more than the shell made of boron/epoxy which Is 
stiffen, but experiences small degree of nonlinearity. 

Figure 4 shows nonlinearity exhibited by the graphite/epoxy cross- 
ply and angle-ply shells under uniform load. Clearly, the angle-ply 
shows greater displacement and also nonlinearity than the cross-ply for 
the same load. 

Figure 5 shows the material behaviour for the boron/epoxy cross-ply 
shell under concentrated load. 

4. Clamped Cylindrical Cross-Ply (0/90) Shell Linder Uniform Load 
The geometry of the shell is shown in Fig. 6. The shell Is madt of 
grpahite/epoxy and the plot of displacement versus load are shown in 
Fig. 6. 


CONCLUSIONS 

A finite element model based on Sander's shell theory, accounting 
for the transverse shear strains is used for the elastic-plastic 
analysis of lamlanted composite shells. The parameters of anisotropy 
reflect the plastic material response by correcting the stress 
components in the Hill's yield function. Numerical results are 
presented for isotropic and laminated shell of cylindrical and >pher1cal 
geometry to demonstrate the validity and efficiency of the present 


Center deflection, w (in.) 


\ Figure 3. Load-deflection curves for a simply supported 

\ spherical shell under uniform transverse load 

[ 



9Q 


100 -I 

— Graph te-epoxv (0°/90°) 


j Graphi te-epoxy (-45°/45‘) 

80 H 



Center deflection, w (in.) 


Figure 4. Load-deflection curves for a simply supported 
spherical shell (see Figure 3 for the geometry) 
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Figure 5. Load-deflection curves for a simply supported spherical shell 
under point load at the center (see Figure 3 for the geometry) 



igure 6. Load-deflection curve for a clamped cy 
for the geometry with L = 50 in., R = 
under uniform transverse load, 
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approach. For the Isotropic case, the present results are In good 
agreement with those available In the literature. 
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APPENDIX I 

Transformation of the stress-strain matrix In Equation (6) 

Let the elastic matrix in the material axes (1,2) be [ D e P ] and '• n 
the body axes be [D ep ] xy 



C 11 

c 12 

C 13 

C 14 

C 15 

[D ep ]j 2 = [Cl = 


°22 

^3 

c 33 

C 24 

C 34 

v 

C 25 

C 35 


sym 



C 44 

C 45 


— 




c 55_ 


"Qii 

Ql2 

q 13 


q i5 



Q 22 

^ q 23 


025 

CD ep ] xy -rqj - 



q 33 

^34 

035 


sym 



Q 44 

q 45 






q 55_ 


then the transformation [4] Is given as* (with m = cose, n = sine) 

Qll = m 4 C 11 + 2m 2 n 2 (C 12 + 2C 33 ) - 4mn(m 2 C 13 + n 2 C 23 ) + n 4 C 22 
Q i 2 = rn 2 n 2 (C 11 + C 22 - 4C 33 ) + 2mn(m 2 - n 2 )(C 13 - C 23 ) + (m 4 + n 4 )C 12 
Q ]3 = m 2 (m 2 - 3n 2 )C 13 + mn[m 2 C 11 - n 2 C„ 2 - (m 2 - n 2 )(C l2 + 2C 66 )} 

+ n 2 (m Z - n 2 )C 26 

^14 ” m ^14 ” mn t ( ^^34 ” ~ ^24 ~ 2 ^ 35 ^ n ^ + ^ 25 n 

Ql5 ~ ^ - mn[(C^ + 2C 3 g)m - (C 2 g + 2C 3 ^)n] - C 2 ^n 

Q 22 = n 4 C u + 2m 2 n 2 (C 12 + 2C 33 ) + 4mn(m 2 C 23 + n 2 C 13 ) + m 4 C 22 
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0 ?3 = m 2 (m 2 - 3 n 2 )C 23 + mn{rrC,j - m 2 C 22 + (m 2 - n 2 ) (C 12 + 2C 33^ 

+ n Z ( 3 m 2 - n 2 )C n 

Q^4 = ^ ^24 mr| [ (^25 + 2 ^34^ m + ^14 + 2 ^35^^ ^15*"* 

^25 “ m ^25 ” mn KC24 _ 2C 3 g)m - (C^g - 2C 3 ^)nJ - C^n 

Q 33 = m 2 n 2 (C 11 + C 22 - 2C 12 ) - 2mn(m 2 - n 2 )(C„ 2 - C 33 ) + (m 2 - n 2 ) 2 C 33 

Q34 = (^34 + (rn - n ) + m n(C^ - ^24^ + ^15 ~ ^ 25 ^ 

Q35 = ^^35 “ ^^34) ~ n ) + m n(C^g ~ "25) + mn ^24 ~ ^ 14 ^ 

Q44 = m?C 44 + 2rnnC 45 + "^55 
^45 * (m 2 - q2 ) C 45 “ mn ( C 44 “ l 55^ 

Q55 = nfl ^ C 55 ~ 2mnC 45 + "^44 

The underscored terms are due to material nonlinearity for an 
orthotropic material. Also note that the constitutive matrix is no 
longer orthotropic. 
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A PPENDIX II 

Stiffness Coefficients 

[K U | = A U [S U | * A 16 ((S 12 1 i- [S 21 !) + Agg[S 22 ] 

+ C 0 (B 16 (S 12 1 + [S 21 ]) + 2B 66 IS 22 ] + C^jglS 22 !) + ^ IS 00 ) 

• J- (« 15 ([S 10 !+[S 01 1 )+A 56 (|S 20 I+(S 02 |)-C o B 56 ([S 20 ]i-(S 02 1 ) 

[K I2 | - A 12 (S 12 | + A 16 [S U I + a 26 is 22 i + a 66 [s 21 1 

+ C o'' 8 26' s22 1 - B 16 IsI1 > - C o 0 66 [s21 D |S °°I 

- < A 14‘ s101 + A 46' s2 °> + C o B 46' s2 °^ 

- ^ ( fl 25 [S ° 2 l + A 56 [s011 - C o B 56 1s01 I> 

[K 13 ] - t- (AjjlS 10 ] + A 16 IS 20 |) * i- (A 12 [S 10 1 + A 26 [S 20 |) + 

Vsf ls2D| + sf [s20 ' - S7 (V s ° 2 i + V s01| > 

- A 14 [S 12 ) + A^iS 11 ] - i£i) ♦ A 46 (S 22 1 

^ 

^ A 56 [S 21 1 , C q[B 46 [S 22 1 . 8 56 [ s21 ]) - A ?5 ig 
IK 14 ] = B n [S U ] + 8 16 ([S 12 ] + [S 21 j) + B,. g [S 22 ] + 

W s21 l + 0 66 l s22 D * itf A 55 [S ‘' 0 l 

* a 15 [s 10 i 4 . A 56 [S 2 °] -^(B 15 [S 01 ) + B 5 .(S 02 1) 

- C 0 B 56 iS 2[ >] 
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|K 15 ] . B, ; ^ 12 ] + 8 16 [S U + B 26 !S Z2 | + B 66 (S 21 | 

t 0 66 lsJ1 » -^ A 45 !S ° 01 
+ A 14 [s10 ' + A 46 IS2 °' ' W{ <B 25 (S 02 1 * 8 56 (S 01 |> 

* 

{k 21 1 = lK 12 j T 

[K 2Z I = A 22 [S 22 ) + A 26 ([S 12 } + is 21 ]) + A g 6 (S U ] - 2 C 0 B 56 [ 5 n ] 

- c 0 ( p 2 6 (?s125 + {s2l]) - c o°66 i5lIi) - ~r ls ° 0) 

“ r 2 ( A 24^ sZ ° 5 + t s ° 2 l) + A 46 US 10 ] * (S 01 ) + C o [S 10 i) 

- C o B 46' s1 °H 

IK 23 } - (A 12 [S 2Q 1 + A 1 Q {S 10 ]) + i- (A 22 [S 2 °] + A 26 [S 10 )} 

' C o^ + lf } (Sl0] -k^ 02 ^^ 0 ^ 

+ A 24 (lS t2 ) - -^2 ( S 00 1 ) + A 25 (S *J + A 4 gt si2 l 
22 

(k 24 i - e 12 {$ 21 ! + s 26 [s 22 ] + 8 16 ts n j + b 66 [s 12 i 

- C o< D 16^ 1 ^ + °66^ 12 ^ ' A 45* S * 

- a 25 is 20 i - A 56 (S 10 ) - <B 14 [S 01 ] ^ a 46 is 02 ]) 
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- W 5 ' 0 ' 


(K 25 3 = B 22 (S 22 ) + B 26 ((S 21 ] + [s 12 l) + rf 66 (S u ! 
- C 0 (D 26 [S 12 J -h 0 66 IS X1 ]) - 1 A 44 fS 00 ! 


* A 24 t s2 °l * A 46 1S10) * SJ ( B 24 ! : ° 2 I * B 46 |s01|) 


C o B 46 1s10 ' 


[K 31 ] - [K 13 1 T 


IK 32 ] = [K 23 ] T 


|K 33 I = A 45 [S 12 1 * A 55 [S U ] * A 44 [S 22 1 * A 45 [S 21 1 


, ,,oo,,l , A 11 , A 12, , 1 Si . A 22„ 

[S »(r“ (p + R ) + R- ' R + R )) 

K 1 K J k 2 k 2 k 1 k 2 


+ ([S 02 l + IS 20 ]) + ^ ( [S 01 ] + [S 10 ]) 

21 21 

+ ^ ((S 02 ] + [S 20 j) + ^ ([S 01 j + [S 10 ]) 
22 22 

+ ^((s 02 ! ♦ [S 20 ]) + ^( |S 01| ♦ [S 10 ]) 


[K 34 ] = A 5e (S 10 ) + A 4g [S 20 ] 

M^^)IS 01 ]M^ + ^)IS 02 I 

K 1 k 2 k 1 k 2 

- + sf> ls °°> + e i4 [s21 ' + 8 ^> s22 ' 

+ B lg (S 11 ] 4 B 56 (S 12 ] 


46 l 
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[K 


35 , 


t . 


i - 


' A 45 [s10) + V s201 

+ ( !li + ! 26 )(S C 1 | + ,!u + !n )ls o 2l 

K 1 K 2 k 1 k 2 



A. . A_. 

, / 14 24 v f r 00 i 

+ ^ _ — ) ( s } + 
K 1 K 2 

b 24 [$ 22 ! + 

b 4 S (s 21 ] 


* b 25 [ S '- ? i . 

► 0 56 ts xl i 



[K 41 ] 

» (K 14 ] T 




[K 42 ] 

- [K 24 ] 1 




[k 43 1 

. (k 34 ] t 




EK 44 ] 

= D 11 IS ii ] + 

0 I 6 ([S 12 ! 

1 + (S 2i |) < 

■ “eel 5 ” 1 + 


+ b 15 ([s 10 ] 

+ E s° 1 1 ) 

+ B 56 ((S 2 °] 

+ ES 02 ]) 

[K 45 ] 

= a 12 f s 12 ] + 

Oi 6 iS ll l 

+ o 26 es 22 i 

+ D 66 [S2!| ' 


+ Bi 4 [S 10 ] 

- B 46 |s2 ° 

1 + 8 25 ts 02 ] 

- b 56 is 01 i 

[K S1 t 

- [ K 15 ] T 




[K 52 ] 

- [ K 25 ! r 




[k 53 I 

- I K 35 ] T 




(k 54 1 

= (K 45 I T 




IK 55 ] 

■ D 22 1s22 > * 

d 26 (Is 12 i 

+ IS 21 ]) + 

° 66 [SU| + 


+ B 24 (IS 20 | 

+ IS 02 ]) 

* B 46 « S1 ° 

1 + [S 01 ]) 


. 00 , 


' 55 ' 


-00 , 


45 


• 00 , 


44 


100 


where 




3ib av. 

^ dx.dx, 

ax ax 1 2 

a s 


-oo 

>ij 


■v^jdx^dx^ 


and the underscored 


terms are due to material nonlinearity. 
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